    pinks = {'#860454','#c51b7c','#dc75ab','#f0b7da','#FF00FF','#ff81c0'};
    greens = {'#276418','#4d921e','#7fbc42','#b9e084','#15b01a','#77AC30'};
    purples = {'#410149','#762a84','#9b6fac','#c1a5cd','#800080','#7e1e9c'};
    browns = {'#523107','#7f3c0a','#bf812c','#e2c17e','#a0522d','#A2142F'};
    pools = {'#003d2e','#01665e','#379692','#81cdc1','#13eac9','#04d8b2'};
    reds = {'#68011d','#b5172f','#d75f4e','#f7a580','#ff0000','#e50000'};
    blues = {'#062e61','#2265ad','#4295c1','#93c5dc','#40e0d0','#00FFFF'};
    
%% Linear polarization
Files_linear = {'data_20230331_100VR_30G_10kHz_REMPI_scan_4', ... % R0, with light shift
   'data_20230509_100VR_30G_10kHz_REMPI_scan_9'}; % R0, no light shift


figure(89)
clf

subplot(1,2,1)
hold on
peak = zeros(1,length(Files_linear));
e_peak = zeros(1,length(Files_linear));
lgs = {'With light shift','No light shift'}
for i = 1:length(Files_linear)
    d = load(Files_linear{i});
    x = d.x;
    ys = d.ys;
    eys = d.es;
    yb = d.yb;
    eyb = d.eb;
    errorbar(x,abs(ys-yb),sqrt(eys.^2+eyb.^2),'o','markersize',8,'linewidth',2,'color', hex2RGB(greens{2*i-1}),'markeredgecolor',hex2RGB(greens{2*i-1}),'markerfacecolor',hex2RGB(greens{4}),'DisplayName',lgs{i})
%     if i == 1 | i==3 | i==7
%          errorbar(x,ys-yb,sqrt(eys.^2+eyb.^2),'-o','markersize',8,'linewidth',2,'color', hex2RGB(greens{1}),'markeredgecolor',hex2RGB(greens{1}),'markerfacecolor',hex2RGB(greens{4}))
%     elseif i == 4 | i==6
%         errorbar(x,ys-yb,sqrt(eys.^2+eyb.^2),'-o','markersize',8,'linewidth',2,'color', hex2RGB(greens{2}),'markeredgecolor',hex2RGB(greens{2}),'markerfacecolor','white')
%     elseif i == 5
%         errorbar(x,ys-yb,sqrt(eys.^2+eyb.^2),'-o','markersize',8,'linewidth',2,'color', hex2RGB(greens{3}),'markeredgecolor',hex2RGB(greens{3}),'markerfacecolor','white')
%     elseif i == 2
%         errorbar(x,ys-yb,sqrt(eys.^2+eyb.^2),'-o','markersize',8,'linewidth',2,'color', hex2RGB(greens{5}),'markeredgecolor',hex2RGB(greens{5}),'markerfacecolor','white')
%     end
    % fitting
%     xft = linspace(min(x),max(x),100);
%     c0 = [100,50,mean(x),10];
%     [c,cerr,fitval] = fit_wrapper(@(b,x)(b(1)*b(2).*2./((4*(x-b(3)).^2+b(2).^2)) + b(4)),c0,x,ys-yb,sqrt(eys.^2+eyb.^2),xft);
%     plot(xft,fitval,'-','linewidth',1.5,'color', hex2RGB(greens{1}))
%     
%     peak(i) = 2*c(1)/c(2)+c(4);
%     e_peak(i) = sqrt(4*c(1)^2/c(2)^2*(cerr(1)^2/c(1)^2+cerr(2)^2/c(2)^2)+cerr(4)^2);
    
        xft = linspace(min(x),max(x),100);
    c0 = [100,50,mean(x)];
    [c,cerr,fitval] = fit_wrapper(@(b,x)(b(1)*b(2).*2./((4*(x-b(3)).^2+b(2).^2)) ),c0,x,abs(ys-yb),sqrt(eys.^2+eyb.^2),xft);
    plot(xft,fitval,'-','linewidth',1.5,'color', hex2RGB(greens{2*i-1}))
    
    peak(i) = 2*c(1)/c(2);
    e_peak(i) = sqrt(4*c(1)^2/c(2)^2*(cerr(1)^2/c(1)^2+cerr(2)^2/c(2)^2));
end

ylim([0 40])
box on; %grid on
xlabel('666 nm laser frequency 449850 GHz + x (MHz)')
ylabel('KRb ion counts')
set(gca,'fontsize',16,'linewidth',1)
set(gcf,'color','white')
legend


if 1
    Files ={'data_20230109_77to110_11.mat' ...
        };

    figure(89)
    
    
    titles = {'KRb N = 0'};
    for i = 1:length(Files)
        load(Files{i});
        d = p.Rb;
        x = d.x_unique;   
        x = x-4000;
    %     y = d.analysis.atomN.val;
    %     err = d.analysis.atomN.err;
        y = d.analysis.peakOD.val;
        err = d.analysis.peakOD.err;
        coeff = max(y);
        y=y/coeff;
        err = err/coeff;
    %     subplot(3,1,i)
        subplot(1,2,2)
        hold on
        errorbar(x,y,err,'o','linewidth',2,'markersize',8,'Displayname',titles{i})
    %     title(titles{i})
    
        
        xft = linspace(min(x),max(x),1000);
        c0 = [-0.5,40,-3580,  -0.5,40,-3540,  -0.5,40,-3500];
        lb = [-20,10,-3600,  -20,10,-3580,  -20,10,-3540];
        ub = [0,200,-3540,  0,200,-3500, 0,200,-3450];
        [c,cerr,fitval] = fit_wrapper2(@(b,x)(1+b(1)*b(2).*2./((4*(x-b(3)).^2+b(2).^2))+b(4)*b(5).*2./((4*(x-b(6)).^2+b(5).^2))+b(7)*b(8).*2./((4*(x-b(9)).^2+b(8).^2)) ...
            ),c0,x,y,err,xft,lb,ub);
        plot(xft,fitval,'-','linewidth',1.5,'displayname','fitting')

        xlabel('666 nm laser frequency 449850 GHz + x (MHz)')
        ylabel('Normalized KRb molecule number')
        box on; %grid on
%         ylim([0 1.05])
%         xlim([-11500 2500])
        set(gca,'fontsize',16,'linewidth',1)
        set(gcf,'color','white')
        legend
    end
end


function [c,cerr,fit_val] = fit_wrapper2(func_hand,c0,x,y,ye,xft,lb,ub)
    if nargin < 5
        ye = 0;
    end
    if nargin < 6
        xft = x;
    end
    opts = optimset('Display','off');
    % [c,~,R,~,~,~,J] = lsqcurvefit(@gaussian2d,c0,xy,im(:),[],[],opts);
    if ye == 0
        [c,R,J] = nlinfit(x,y,func_hand,c0,opts); % slightly faster
    else
%         [c,R,J] = nlinfit(x,y,func_hand,c0,opts,'weights',1./ye.^2);
         [c,~,R,~,~,~,J] = lsqcurvefit(func_hand,c0,x,y,lb,ub);
    end
    ci = nlparci(c, R,'jacobian',J,'alpha',.3173); % 1 sigma
    cerr = abs(ci(:,2)-ci(:,1))'/2;
    fit_val = func_hand(c,xft);
end


function color = hex2RGB(hex)
    color = sscanf(hex(2:end),'%2x%2x%2x',[1 3])/255;
end

function [c,cerr,fit_val] = fit_wrapper(func_hand,c0,x,y,ye,xft)
    if nargin < 5
        ye = 0;
    end
    if nargin < 6
        xft = x;
    end
    opts = optimset('Display','off');
    % [c,~,R,~,~,~,J] = lsqcurvefit(@gaussian2d,c0,xy,im(:),[],[],opts);
    if ye == 0
        [c,R,J] = nlinfit(x,y,func_hand,c0,opts); % slightly faster
    else
        [c,R,J] = nlinfit(x,y,func_hand,c0,opts,'weights',1./ye.^2);
    end
    ci = nlparci(c, R,'jacobian',J,'alpha',.3173); % 1 sigma
    cerr = abs(ci(:,2)-ci(:,1))'/2;
    fit_val = func_hand(c,xft);
end